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This  paper  provides  references  for  the  Navy’s  existing 
models  for  shallow-water  propagation.  The  shallow-water  acoustic 
environment  is  then  briefly  described,  followed  by  a description 
of  sound  propagation.  Four  basic  models  of  sound  propagation  in 
the  ocean  most  relevant  to  shallow-water  propagation  are  derived 
from  their  common  wave  equation  origin:  ray  theory,  spectral 
theory,  normal  mode  theory,  and  the  parabolic  equation  method. 
Some  results  from  these  models  are  discussed. 

INTRODUCTION 

The  Oceanographic  and  Atmospheric  Master  Library1 
(OAML)  contains  a description  of  Navy  models  and 
databases.  The  Navy’s  use  of  shallow-water  models  can  be 
summarized  as  follows:  Presently,  the  Navy’s  standard 
shallow- water , passive  acoustic  propagation  model,  Collosus 
II,  is  an  empirical  fit1  to  observed  data.  It  is  used  for  areas 
where  no  significant  environmental  data  are  available; 
otherwise  RAYMODE1  is  used  for  range-independent 
environments.2 

For  range-dependent  environments,  the  parabolic 
equation  (PE)  models3  are  used,  including  the  standard  split- 
step  and  the  more  recent  finite-element  (FEPE)  model. 
ASTRAL1  is  used  for  higher  frequencies,  say  above  250 
Hz.  The  biggest  problem  facing  shallow-water  propagation 
prediction  is  the  adequacy/resolution  of  environmental 
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databases.  Various  compilations  of  existing  models  are 
found  in  Refs.  1,  4,  and  5.  The  general  subject  of  ocean 
acoustic  modeling  and  numerical  methods  is  covered  in  Ref. 
3.  The  purpose  of  this  article  is  to  briefly  review  the 
existing  shallow-water  computer  models  from  the  point  of 
view  of  their  physics  and  their  relationship  to  each  other. 
The  frequency  regime  covered  here  is  up  to  a few  kilohertz. 
High  frequencies  do  not  propagate  very  long  distances,  and 
hence  are  not  dependent  on  gross  shallow- water  geophysical 
properties  that  comprise  the  acoustic  waveguide.  Reference 
6 is  a compilation  of  many  aspects  of  high-frequency 
acoustics. 

Although  shallow-water  acoustics  has  received  great 
attention  recently,  there  already  is  a rich  baseline  knowledge 
of  acoustic  propagation.  Great  progress  has  been  made  in 
understanding  this  subject  since  die  Bradley-Urick7  compila- 
tion of  transmission  loss  for  various  shallow-water  scenarios 
appeared,  The  plots  in  this  latter  reference  showed  a 40  dB 
spread  at  10  Ian  and  a 50  dB  spread  at  100  km.  In  this 
paper,  we  review  the  progress  in  modeling  sound  propaga- 
tion in  shallow  water  and  the  subsequent  understanding  that 
modeling  and  data  analysis  has  provided.  The  bottom  line 
is  that  the  environmental  inputs  to  acoustic  models  remain 
the  crucial  factor  for  quantitatively  predicting  propagation 
conditions. 

We  start  this  paper  with  a qualitative  description  of 
sound  propagation  in  shallow  water.  Next  we  describe  the 
generic  propagation  models  in  enough  detail  to  give  the 
reader  a sense  of  both  their  common  underlying  physics  and 
their  different  theoretical  and  algorithmic  approaches.  We 
conclude  the  paper  with  illustrative  results  from  the  various 
models. 
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DESCRIPTIVE  SHALLOW-WATER  ACOUSTICS 

In  this  section  we  describe  the  basic  features  of 
shallow-water  environments  and  how  these  features  control 
sound  propagation.3*8  9 See  the  Appendix  for  discussion  of 
units. 

Environment 

By  shallow  water,  this  author  means  continental  shelf 
and  slope  such  that  sound  propagation  is  best  described  as 
a waveguide  phenomenon.  The  waveguide  is  defined  by  the 
upper  air-sea  boundary  and  the  ocean  bottom  interface, 
which  in  a more  general  sense  may  be  a composite  of 
stratified  and/or  nonstratified  layers.  The  latter  is  often 
called  a range-dependent  environment  that  also  includes  the 
case  of  the  sound  speed  varying  as  a function  of  horizontal 
position,  not  just  depth.  The  acoustic  modeling  significance 
of  a range-dependent  environment  is  that  the  wave  equation 
describing  acoustic  propagation  is  not  a " separable"  partial 
differential  equation,  and  therefore  usually  requires  more 
computationally  intensive  methods  than  separable  equations. 

The  two  most  common  classes  of  sound  speed  profiles 
in  shallow  water  are  schematically  represented  in  Fig.  1. 
Typically,  much  of  what  is  classified  as  shallow  water  is  not 
deep  enough  for  the  depth-pressure  term  in  the  temperature- 
salinity-pressure  dependence  of  sound  speed  to  be  signifi- 
cant. Thus  the  winter  profile  tends  to  isovelocity  simply 
because  of  mixing,  whereas  the  summer  profile  has  a higher 
sound  speed  near  the  surface  due  to  heating.  The  sound 
speed  profile  may  vary  spatially  and  temporally.  Often,  the 
temporal  variations  are  caused  by  internal  wave  fields  that 
tend  to  be  more  ordered  in  shallow  water  than  in  deep 
water.  Sometimes  confined  internal  wave  packets,  referred 
to  as  solitons,  produce  localized,  travelling  ocean  distur- 
bances. 

The  ocean  bottom  has  a great  impact  on  sound  propa- 
gation due  to  its  material  properties.  Table  1 summarizes 
some  of  the  acoustic  properties  of  the  ocean  bottom.3  The 
effects  of  these  bottom  properties  are  briefly  described  in 
the  next  subsection. 

Acoustics 

Much  of  the  qualitative  behavior  of  sound  in  shallow 
water  can  be  described  by  a combination  of  a simple 
application  of  Snell's  law  and  some  understanding  of  the 
reflection  properties  of  the  ocean  bottom.  Basically,  Snell's 
law  states  that  sound  bends  toward  regions  of  low  sound 
speed.  Hence,  as  described  above,  summer  sound  speed 


profiles  produce  rays  that  bend  more  toward  the  bottom 
than  winter  profiles  in  which  the  rays  tend  to  be  straight. 
This  implies  two  effects  with  respect  to  the  ocean  bottom: 

• for  a given  range,  there  are  more  bounces  off  the 
ocean  bottom  in  summer  than  in  the  winter;  and 

• the  ray  angles  intercepting  the  bottom  are  steeper  in 
the  summer  than  in  the  winter. 

A qualitative  understanding  of  the  reflection  properties  of 
the  ocean  bottom  should  therefore  be  very  revealing  of 
sound  propagation  in  summer  vs  winter.  Basically,  near- 
grazing incidence  is  much  less  lossy  than  larger,  more 
vertical  angles  of  incidence.  Because  summer  condition 
propagation  paths  have  more  bounces,  each  of  which  are  at 
steeper  angles  than  winter  paths,  we  can  expect  summer 
shallow-water  propagation  to  be  much  worse  than  in  winter 
conditions,  which  indeed  they  are.  This  result  is  often 
tempered  by  rough  winter  sea  surface  conditions  that 
generate  large  scattering  losses  at  the  higher  frequencies. 

Bottom  Loss 

Modeling  the  shallow-water  environment  as  a wave- 
guide yields  a qualitative  description  of  sound  propagation 
that  is  very  helpful  in  understanding  more  sophisticated 
model  outputs.  However,  to  understand  the  effects  of  the 
ocean  bottom,  let  us  first  review  the  reflection  properties  of 
an  interface  separating  two  media. 

Ocean  bottom  sediments  are  often  modeled  as  fluids 
since  the  rigidity  (and  hence  the  shear  speed)  of  the  sedi- 
ment is  usually  considerably  less  than  that  of  a solid,  such 
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Table  1 — Geoacoustic  Properties  of  Shallow  Water  and  Continental  Slope  Ocean  Bottoms 
(fe  is  the  depth  below  the  water-bottom  interface) 


Bottom 

Type 

p 

(%) 

pb/pw 

(m/s) 

C. s 

(m/s) 

(dB/X,) 

a, 

(d  B/X,) 

Clay 

70 

1.5 

1.00 

1500 

<100 

0.2 

1.0 

Silt 

55 

1.7 

1.05 

1575 

«(1) 

1.0 

1.5 

Sand 

45 

1.9 

1.1 

1650 

P (2) 

0.8 

2.5 

Gravel 

35 

2.0 

1.2 

1800 

r (3) 

*"  i 

0.6 

1.5 

Moraine 

25 

2.1 

1.3 

1950 

600 

0.4 

1.0 

Chalk 

— 

2.2 

1.6 

2400 

1000 

0.2 

0.5 

Limestone 

— 

2.4 

2.0 

3000 

1500 

0.1 

0.2 

Basalt 

— 

2.7 

3.5 

5250 

2500 

0.1 

0.2 

c,(l)  = 80  z°'3  c/3)  = 180  z°'3 

cs<2)  = 100  z°- 3 cw  = 1500  m/s,  pK  = 1000  kg/m3 


as  rock.  In  the  latter  case,  which  applies  to  the  "ocean 
basement"  or  the  case  where  there  is  no  sediment  overlying 
the  basement,  the  medium  must  be  modeled  as  an  elastic 
solid.  This  means  that  it  supports  both  compressional  and 
shear  waves. 

Reflectivity,  the  amplitude  ratio  of  reflected  and 
incident  plane  waves  at  an  interface  separating  two  media, 
is  an  important  measure  of  the  effect  of  the  bottom  on 
sound  propagation.  For  an  interface  between  two  fluid  semi- 
infinite  halfspaces  with  density  p,  and  sound  speed  i = 
1,2,  and  a unit  amplitude  incident  wave  of  the  form  exp 
[i(k1  r + kh  z)  - /cor]  [Fig.  2(a)],  the  reflectivity  is  given 
by 

m = — — — ♦ (i) 

pK  + p\K 

with 

ku  = (co/c.)  sin  6i  ■ kt  sin  0(.;  i = 1,2.  (2) 


The  incident  and  transmitted  grazing  angles  are  related  by 
Snell’s  law, 

k±  = kx  cos  01  = k2  cos  02 , (3) 

where  the  incident  grazing  angle  6X  is  also  equal  to  the  angle 
of  the  reflected  plane  wave.  R(B)  is  also  referred  to  as  the 
Rayleigh  reflection  coefficient,  and  has  unit  magnitude  (total  \ 
internal  reflection)  when  the  numerator  and  denominator  of 
Eq.  (1)  are  complex  conjugates.  This  occurs  when  k ^ is 
purely  imaginary.  Using  Snell’s  law,  Eq.  (3),  for  determin- 
ing S2  in  terms  of  the  incident  grazing  angle,  we  obtain  the 
critical  grazing  angle  below  which  there  is  perfect 
reflection, 

cos  6c  = c/c2,  (4) 

so  that  a critical  angle  can  exist  only  when  the  speed  in  the 
second  medium  is  higher  than  that  of  the  first.  Below  the 
critical  angle,  the  reflected  wave  undergoes  a change  of 
phase  e,  given  by  the  argument  of  the  complex  reflection 
coefficient;  in  this  region  [and  using  Eq.  (2)],  the  reflection 
coefficient  can  be  written  as  a complex  number  with  unit 
magnitude 
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R(0)  = 


where 


e - 2 tan' 


,/cos2  0, 

- (c2Jcl) 

/<; 

cos2  0, 

./cos2  0, 

- (c2/c2) 

/l  - 

cos2  0, 

, P,/C°s2 

0,  - ichcb 

Pi  / 

1 - cos2  0, 

= e",  (5) 


(6) 


Note  that  p2  = 0 is  the  pressure  release  case,  which  is  a 
good  approximation  to  that  water/air  interface;  for  this  case 
there  is  a 180-deg  phase  change.  For  the  two-fluid  interface 
under  consideration,  the  phase  change  goes  from  180  deg  at 
horizontal  grazing  incidence  to  no  phase  change  at  the 
critical  angle.  On  the  other  hand,  when  the  density  and 
speed  of  the  second  medium  is  very  large,  corresponding  to 
the  rigid  interface  condition,  there  is  no  phase  change.  Since 
it  turns  out  that  most  shallow-water  propagation  corresponds 
to  paths  with  very  low  grazing  angles,  the  pressure  release 
condition  is  actually  a better  approximation  of  a two-fluid 
interface  than  the  rigid  boundary  condition. 


b) 


Fig.  2 — The  reflection  and  transmission  process 


with  the  total  impedance  of  the  second  medium  being 


Z*  - Z*  sin2  2 e„  + cos2  20*.  (9) 


Using  Eq.  (2),  Eq.  (1)  can  also  be  rewritten  as 


Snell’s  law  for  this  case  is 


R(6)  = — Cz/Sin  e2  ~ P|  C|/Sin  6-l  m h ~ Z|  (7) 
p2  cj sin  02  + px  c/sin  01  + Zx 

where  the  right-hand  side  is  in  the  form  of  impedances  Z, 
(0.)  = p.c/sin  0f,  which  are  the  ratios  of  the  pressure  to  the 
vertical  particle  velocity  at  the  interface  in  the  ilh  medium. 
Written  in  this  form,  more  complicated  reflection  coeffi- 
cients become  intuitively  plausible.  Consider  the  case  in 
Fig.  2(a)  in  which  the  second  medium  is  elastic  and  thus 
supports  shear  as  well  as  compressional  waves  with  sounds 
speeds  and  c2p,  respectively.  The  Rayleigh  reflection 
coefficient  is  then  given  by 

R(6)  = Z*L Z‘ , (8) 

Z*.  + Z, 


kx  cos  0,  = cos  0^  = kty  cos  0^.  (10) 

Shear  speeds  of  less  than  a few  hundred  m/s  have  very 
little  effect  on  sound  propagation.  However,  sediments  with 
shear  speeds  greater  than  a few  hundred  m/s  but  still  much 
lower  than  the  speed  of  sound  in  water  (say,  up  to  1000 
m/s)  appear  to  be  highly  absorbing  because  the  propagation 
of  shear  waves  is  another  degree  of  freedom  for  sound  to  be 
transmitted  away  from  the  water  column.  For  even  greater 
shear  speeds  but  still  less  than  the  speed  of  sound  in  water, 
the  bottom  appears  almost  transparent  and  waveguide 
propagation  is  severely  attenuated.  For  the  case  of  the  shear 
speed  exceeding  the  sound  speed  in  water,  the  shear  speed 
is  the  dominant  bottom  sound  speed  parameter  because  its 
value  is  closest  to  the  water  sound  speed  (the  bottom 
compressional  speed  must  be  larger  than  the  shear  speed). 
For  this  case,  a reasonable  approximation  (which  also 
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neglects  the  existence  of  interface  waves)  to  the  magnitude 
of  the  bottom  reflectivity  is  to  treat  the  ocean  bottom  as  an 
effective  fluid  whose  compressional  speed  is  the  shear  speed 
of  the  elastic  bottom.  Thus,  the  ratio  of  the  shear  speed  to 
the  water  sound  speed  will  determine  the  critical  angle. 

In  lossy  media,  attenuation  can  be  included  in  the 
reflectivity  formula  by  taking  the  sound  speed  as  complex 
so  that  the  wavenumbers  are  subsequently  also  complex, 
-*k(  • f a*.  Figure  2(b)  depicts  a simple  bottom  loss  curve 
derived  from  the  Rayleigh  reflection  coefficient  formula. 
Both  the  densities  and  sound  speed  of  the  second  medium 
j are  larger  than  those  in  the  first  medium,  with  unit 
reflectivity  indicating  perfect  reflection.  The  phase  change 
upon  reflection  is  also  shown.  For  loss  in  dB,  0 dB  is 
perfect  reflecting,  6 dB  loss  is  an  amplitude  factor  of  one- 
half,  12  dB  loss  is  one-fourth,  etc.  For  a lossless  bottom, 
severe  loss  occurs  above  the  critical  angle  in  the  water 
column  due  to  transmission  into  the  bottom.  For  the  lossy 
(more  realistic)  bottoms,  only  partial  reflection  occurs  at  all 
angles.  With  paths  involving  many  bottom  bounces 
(shallow-water  propagation),  bottom  losses  as  small  as  a 
few-tenths  of  a dB  per  bounce  accumulate  and  become 
significant  because  the  propagation  path  may  involve  many 
tens  of  bounces.  A description  of  a reflectivity  database  is 
given  in  Section  2.2.8  of  Ref.  1, 

Qualitative  Description  of  Waveguide  Propagation 


SOUND 

SPEED 

PROFILE 


RANGE 


OCEAN  BOTTOM 


(a)  Long-distance  propagation  occurs  within  a cone  of  2 0C 


b)  RAY-MODE  ANALOGY 


(b)  Geometry  for  the  constructive  interference  of  wavefronts 
to  form  a mode 

Fig.  3 — Ocean  waveguide  propagation 


For  simplicity,  we  consider  an  isovelocity  waveguide 
bounded  above  by  the  air/water  interface  and  below  by  a 
two-fluid  interface.  Hence,  we  will  have  perfect  reflection, 
with  a 180-deg  phase  change  at  the  surface.  For  paths  more 
horizontal  than  the  bottom  critical  angle,  there  will  also  be 
perfect  reflection  with  an  angle-dependent  phase  change 
given  by  Eq.  (6).  Therefore,  as  schematically  indicated  in 
Fig.  3(a),  ray  paths  within  a cone  of  2 6C  will  propagated 
unattenuated  down  the  waveguide.  Because  the  up  and  down 
going  rays  have  equal  amplitudes,  preferred  angles  will 
exist  such  that  perfect  constructive  interference  can  occur. 
These  particular  angles  can  be  associated  with  the  normal 
modes  of  the  waveguide  as  formally  derived  from  the  wave 
equation  in  the  next  section.  However,  it  is  instructive  to 
understand  the  geometric  origin  of  the  waveguide  modal 
structure. 

Figure  3(b)  is  a schematic  of  a ray  reflected  from  the 
bottom  and  then  the  surface  of  a "Pekeris"  waveguide  (an 


environment  with  constant  sound  speeds  and  densities  in  the 
water  column  and  fluid  bottom,  respectively).  Consider  a 
ray  along  the  path  ACDF  and  its  wavefront,  which  is 
perpendicular  to  the  ray.  The  two  down-going  rays,  AC  and 
DF,  which  are  of  the  form 

<t>dowr,,up  * exp  li(k±  r±kllz)-iut ] (H) 

will  constructively  interfere  if  points  B and  E have  a phase 
difference  of  an  integral  number  of  360  deg  (and  similarly 
for  up-going  rays).  The  acoustic  phase  change  along  BCDE 
is  given  by  the  product  of  the  wavenumber  and  the  distance 
along  the  ray  plus  the  change  of  phase  at  C and  a 180-deg 
phase  change  at  D,  the  pressure  release  surface. 

Noting  that  ~KB  = ZXE,  the  condition  for  constructive 
interference  is 
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k 


d 

sin  eK 


— C- — cos(ir  - 20  ) 
sin  6 . 


+ £ + 7 r = 2(n  - 1)  7r;  n = 1,2,..., 
which,  after  using  Eq.  (6),  reduces  to 


tan 


(12) 


(13) 


where  kn  is  the  horizontal  wavenumber,  kL , corresponding 
to  the  discrete  angle  6„  and  we  have  converted  from  angles 
to  wavenumbers  by  using  the  relations  c,/c2  = k2/kl  and 

*n  = sin  0„  = - k2n  . 


Equation  (13)  is  a transcendental  equation  in  k„  with 
real  roots  in  the  range 

k2  < k < k.  ; (14) 

2 n l 

kn  can  be  thought  of  as  a wavenumber  of  a horizontal  wave 
traveling  down  the  waveguide  with  a phase  velocity  of  cn  m 
o)fkn.  Figure  3(b)  shows  that  the  more  vertical  the  ray,  the 
more  horizontal  the  wavefront.  Hence,  a vertical  ray 
corresponds  to  a horizontal  wavefront,  and  the  horizontal 
phase  velocity  along  such  a wavefront  must  be  infinite  since 
the  wavefront  has  constant  phase.  However,  Ineq.  (14) 
limits  the  phase  velocity  to  a maximum  of  c2,  that  is,  rays 
more  vertical  than  the  critical  angle  do  not  propagate  (more 
than  a few  bounces)  down  the  waveguide.  (For  the  simple 
limiting  case  of  a waveguide  bounded  above  and  below  by 
pressure  release  surfaces,  the  rays  do  approach  the  vertical 
so  that  the  horizontal  phase  velocity  cn  goes  to  infinity  with 
increasing  mode  number  n.)  On  the  other  hand,  a horizontal 
ray  has  a vertical  wavefront,  and  so  the  phase  is  constant  in 
the  vertical.  In  this  case,  the  horizontal  phase  velocity  is 
limited  to  the  medium  speed.  The  vertical  amplitude  of  the 
wave  can  be  obtained  by  looking  at  the  z part  of  the  fields 
described  by  Eq.  (11).  The  field  for  a particular  n will  be 
a superposition  of  the  up  and  downgoing  fields  of  Eq.  (11). 
However,  we  need  a negative  sign  between  the  two 
amplitudes  so  that  the  field  vanishes  at  the  surface.  The 
result  is  that  the  superposition  of  discrete  up  and  downgoing 
waves  results  in  a vertical  amplitude  distribution  in  the 
waveguide  of  the  form. 


4>n  ~ w„(z);  un(z)  = sin  (Jk?  - k2„  z)  . 


(15) 


The  w„’s  are  called  the  normal  modes  of  the  waveguide, 
and  they  each  propagate  with  phase  velocity  c„.  The  total 
field  in  the  waveguide  is  the  sum  of  all  the  normal  mode 
terms  in  Eq.  (15).  The  vertical  distribution  can  be  thought 
of  as  a superposition  of  up  and  downgoing  plane  waves  at 
discrete  propagation  angles  within  the  cone  ±6C . This 
discussion  is  qualitatively  correct  in  general,  but  is 
quantitatively  limited  to  an  isovelocity  water  column 
overlying  a fluid  bottom  with  a constant  soundspeed,  i.e., 
a Pekeris  waveguide.  In  the  next  section  we  show  that  Eq. 
(13)  is  the  eigenvalue  equation  that  comes  from  normal 
mode  theory.  If  the  bottom  medium  is  attenuating,  then  the 
wavenumber  in  the  bottom  will  be  complex.  In  this  case, 
the  individual  eigenvalues  kn  become  complex,  and  there  is 
an  additional  factor  of  exp(-anr)  multiplying  Eq.  (15). 
Physically,  the  normal  mode  attenuation  coefficients  an 
increase  with  n since  the  higher  order  modes  are  more 
vertical  and  therefore  correspond  to  rays  that  have  greater 
loss  per  bounce  and  more  bounces  per  distance  down  the 
waveguide. 

A simple  limiting  case  is  when  we  approximate  the 
bottom  as  pressure  release.  The  density  of  the  second 
medium  vanishes  so  that  the  right-hand  side  of  Eq.  (13)  is 

zero.  This  reduces  to  d\jkf  - k2n  = nir  so  that  the  normal 
mode  functions  of  this  waveguide,  the  7 r phase-shifted  sum 
of  the  up  and  down  going  waves  of  Eq.  (11),  is  of  the  form 
s\n(mrz/d),  which  are  a set  of  orthogonal  functions  that 
vanish  at  both  surfaces. 

SOUND  PROPAGATION  MODELS 

Sound  propagation  in  the  ocean  is  mathematically 
described  by  the  wave  equation,  whose  parameters  and 
boundary  conditions  are  descriptive  of  the  ocean 
environment.3  There  are  essentially  four  types  of  models 
(computer  solutions  to  the  wave  equation)  to  describe  sound 
propagation  in  the  sea:  ray  theory,  the  spectral  or  fast  field 
program  (FFP),  normal  mode  (NM),  and  parabolic  equation 
(PE).  All  of  these  models  allow  for  the  fact  that  the  ocean 
environment  varies  with  depth.  A model  that  also  takes  into 
account  horizontal  variations  in  the  environment  (i.e., 
sloping  bottom  or  spatially  variable  oceanography)  is  termed 
range-dependent.  For  high  frequencies  (a  few  kilohertz  or 
above),  ray  theory  is  the  most  practical.  The  other  three 
model  types  are  more  applicable  and  usable  at  lower 
frequencies  (below  a kilohertz).  The  hierarchy  of 
underwater  acoustics  models  is  shown  in  schematic  form  in 
Fig.  4. 
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the  limiting  case  of  a fluid  medium,  where  shear  waves  do 
not  exist,  §>(r,  z ) represents  the  compressional  potential  0(r, 
z).  Most  underwater  acoustic  applications  involve  only 
compressional  sources,  which  only  excite  the  P and  SV 
potentials,  eliminating  the  SH  potential  A(r,  z).  The 
displacement  potentials  satisfy  the  Helmholtz  equation  with 
the  appropriate  compressional  or  shear  sound  speeds  cp  or 
respectively. 


spectral 

Cn  = 

A + 2 p 

1/2 

> Cs  = 

Fig.  4 — Hierarchy  of  underwater  acoustics  models 
(note:  *'TD*’  refers  to  time  domain) 

p 

P 

_Pm 

1/2 


(19) 


where  X and  p are  the  Lame  constants. 


The  Wave  Equation  and  Boundary  Conditions 

The  wave  equation  is  typically  written  and  solved  in 
terms  of  pressure,  displacement,  or  velocity  potentials.  For 
a velocity  potential  <p9  the  wave  equation  in  cylindrical 
coordinates  with  the  range  coordinates,  denoted  by  r = 
(x,y),  and  the  depth  coordinate  denoted  by  z (taken  positive 
downward)  for  a source  free  region  is 


The  most  common  plane  interface  boundary  conditions 
encountered  in  underwater  acoustics  are  described  below: 
For  the  ocean  surface,  there  is  the  pressure  release 
condition  where  the  pressure  (normal  stress)  vanishes;  for 
the  appropriate  solution  of  the  Helmholtz  equation,  this 
conditions  is 


p ~0,<p  =0  or  0=0. 


(20) 


VV(r ,z,t)  - 


j_  ay(r,z,o 

c2  dt2 


= 0, 


(16) 


where  c is  the  sound  speed  in  the  wave  propagating 
medium.  With  respect  to  the  velocity  potential,  the  velocity 
v and  pressure  p are  given  by 


v = v<p,  p = - p 


dtp 

~dt' 


(17) 


where  p is  the  density  of  the  medium.  The  wave  equation 
is  most  often  solved  in  the  frequency  domain,  that  is,  a 
frequency  dependence  of  exp(-iw/)  is  assumed  to  obtain  the 
Helmholtz  equation  ( K s o >/c), 


V*<p(r9z)  +K2(p(r,z)  = 0. 


(18) 


The  interface  between  the  water  column  (layer  1)  and  an 
ocean  bottom  sediment  (layer  2)  is  often  characterized  as  a 
fluid-fluid  interface.  The  continuity  of  pressure  and  vertical 
particle  velocity  at  the  interface  yields  the  following 
boundary  conditions  in  terms  of  pressure: 


Pi 


or  velocity  potential: 


P,  ~dz 


1 bl± 

P1 


(21) 


P,<Pi  = Pi*v 


3 <Pi 
~dz 


= ^ (22) 

dz  ‘ 


In  underwater  acoustics,  both  fluid  and  elastic  (shear 
supporting)  media  are  of  interest.  In  elastic  media,  the  field 
can  be  expressed  in  terms  of  three  scalar  displacement 
potentials,  0(r,  z),  0 (r,  z),  and  A (r,  z),  corresponding  to 
compressional  (P),  vertically  polarized  (SV),  and 
horizontally  polarized  shear  waves  (SH),  respectively.  In 


These  boundary  conditions  applied  to  the  plane  wave  fields 
in  Fig.  2(a)  yield  the  Rayleigh  reflection  coefficient  given 
by  Eq.  (1). 

For  an  interface  separating  two  solid  layers,  the 
boundary  conditions  are  continuity  of  vertical  displacement 
w(r,  Z/),  tangential  displacements  u(r,  z,)  = ( uxi  uy),  normal 

In 


stress  n ~ and  tangential  stresses  t = (oxl 
cylindrical  coordinates  with  azimuthal  symmetry,  the  radial 


UNCLASSIFIED 


50th  Anniversary  Issue  281 


UNCLASSIFIED 


UNCLASSIFIED 


1062 


UNCLASSIFIED 

WILLIAM  A.  KUPERMAN 


and  vertical  components  of  displacements  (in  homogeneous 
media)  u and  w,  respectively,  are 


u(r,z)  = 


d<t> 

dr 


3V 

dz2’ 


(23) 


Ray  Theory 

Ray  theory10  is  a geometrical,  high-frequency  approx- 
imate solution  to  Eq.  (27)  of  the  form 


G(R)  = A(R)  exp[/S(R)],  <28> 


w(r,z)  - 0L  - i ± r at, 

dz  r dr  dr 


and  the  normal  and  tangential  stresses  are 


t v ~ v dw  , du 

oJr,z)  = (X  + 2\i)  — + X— , 


(24)  where  the  exponential  term  allows  for  rapid  variations  as  a 
function  of  range,  and  ^4(R)  is  a more  slowly  varying 
"envelope"  that  incorporates  both  geometrical  spreading  and 
loss  mechanisms. The  geometrical  approximation  is  that  the 
amplitude  varies  slowly  with  range  (i.e.,  (1  lA)^1  A < < 
j K2),  so  that  Eq.  (27)  yields  the  eikonal  equation 

(25) 

(VS)2  = K2.  <29) 


°r,z  = M 


du  + dw 
dz  dr 


The  ray  trajectories  are  perpendicular  to  surfaces  of  constant 
phase  (wavefronts)  S,  and  may  be  expressed  mathematically 
as 


Continuity  of  these  quantities  at  the  interface  between  two 
solids  are  the  boundary  conditions.  For  a fluid-solid 
interface,  the  rigidity  n vanishes  in  the  fluid  layer,  and  the 
tangential  stress  in  the  solid  layer  vanishes  at  the  boundary. 
If  at  least  one  of  the  media  is  elastic,  these  boundary 
conditions  permit  the  existence  of  interface  or  surface 
waves,  such  as  Rayleigh  waves,  at  the  interface  between  a 
solid  and  vacuum,  Scholte  waves  at  a fluid-solid  interface, 
and  Stoneley  waves  at  a solid-solid  interface.  These  waves 
are  excited  when  the  source  is  acoustically  close,  in  terms 
of  wavelengths,  to  the  interface. 

The  Helmholtz  equation  for  an  acoustic  field  from  a 
point  source  with  angular  frequency  co  is 


V*G(r,z)  + K\r,z)  = - 52(r  - r,)  6(z  - z)\ 


K2(r,z) 


c2(r,z)’ 


(27) 


where  the  subscript  "s"  denotes  the  source  coordinates.  The 
range-dependent  environment  manifests  itself  as  the 
coefficient  K ? (r,z)  of  the  partial  differential  equation  for  the 
appropriate  sound  speed  profile.  The  range-dependent 
bottom  type  and  topography  appear  as  boundary  conditions 
on  scalars  and  tangential  and  normal  quantities,  as  discussed 
above. 


where  / is  the  arc  length  along  the  direction  of  the  ray  and 
R is  the  displacement  vector.  The  direction  of  average  flux 
(energy)  follows  that  of  the  trajectories,  and  the  amplitude 
of  the  field  at  any  point  can  be  obtained  from  the  density  of 
rays.  See  chapter  3 of  Ref.  3 for  more  details. 

The  ray  theory  method  is  computationally  rapid  and 
extends  to  range-dependent  problems.  Furthermore,  the  rays 
give  a physical  picture  of  the  acoustic  paths.  It  is  helpful  in 
describing  how  noise  redistributes  itself  when  propagating 
long  distances  over  paths  that  include  shallow  and  deep 
environments  and/or  mid-latitude  to  polar  regions.  The 
disadvantage  of  conventional  ray  theory  is  that  it  does  not 
include  diffraction  and  such  effects  that  describe  the  low- 
frequency  dependence  ("degree  of  trapping")  of  ducted 
propagation.  In  shallow  water,  many  rays  with  many  bottom 
bounces  are  required.  Hence,  small  errors  from  the  high- 
frequency  approximation  will  accumulate  rapidly  with 
range. 

Full  Field  Solution  or  "Fast  Field  Program  (FFP)’T 

The  wave  equation  can  be  solved  efficiently  with 
spectral  methods3  when  the  ocean  environment  does  not 
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vary  with  range.  The  term  "Fast  Field  Program"  is  used 
because  the  spectral  methods  became  practical  with  the 
advent  of  the  fast  Fourier  transform  (FFT).  Assuming  a 
solution  of  Eq.  (27)  of  the  form 

G<r,z)  = ” d2kg(k,z,z) 

2tt  J-«  (31) 

exp[ik<r  - rs)] 

leads  to  an  equation  for  the  depth-dependent  Green’s 
function  g(k,  z,  zs), 

+ (K2(z)  - k2)  s = - — 8(z  - z.).  (32) 


structure  of  g(k,z,zs)  and  the  normal  modes  of  the 
waveguide.  This  spectral  solution  technique  has  recently 
been  extended  to  include  range-dependent  environments.13 

For  completeness,  we  mention  an  alternative  method14*15 
to  evaluate  Eq.  (31)  using  a ray  like  representation  solution 
[Eq.  (32)].  Up  and  downgoing  solutions  that  satisfy  the 
upper  and  lower  boundary  conditions  are  constructed  by 
using  the  appropriate  boundary  reflection  coefficients.  The 
result  is  an  infinite  sum  of  multiple  bouncing  rays.  Although 
often  such  solutions  are  awkward  for  shallow-water 
problems  (too  many  bounces  combined  with  approximation 
errors),  circumstances  that  involve  reverberation  formulated 
in  terms  of  scattering  at  specific  bounce  locations  can  be 
handled  with  such  a technique. 


Normal  Mode  Model  (NM) 


Furthermore,  we  assume  azimuthal  symmetry  kr  > 2i r and 
rs  = 0 so  that  Eq.  (31)  reduces  to  (after  integrating  over  the 
azimuthal  angle  and  converting  the  Bessel  function  to 
Hankel  functions  and  using  the  asymptotic  form  of  the 

Hankel  function  (HQl  (kr)  -+  \J(2hrkr)  expfifcr  - 7t/4)]) 


Gjr,z)  = exp^— — -4^-  f °° dk(k)mg(k9z,z)ex p(ikr).  (33) 
(2rr)m  J° 


We  now  convert  the  above  integral  to  an  FFT  form  by 
setting  km  = k0  + mAk;  rn  = r0  tiAr  where  n9  m = 0, 
1,..  JV-1  with  the  additional  condition  ArAk  = 2 n/N where 
N is  an  integral  power  of  two: 


Akexp[i(k(jrn  - *74) 

G(v) ' — (2^r — 

A-l 

Xmexp(2irimn/N) 


Xm  = eXP  KM)’ 

Although  the  method  is  labeled  "fast  field,"  it  is  fairly  slow 
because  of  the  time  required  to  calculate  the  Green’s 
functions  [solve  Eq.  (32)]  many  times.  However,  it  has 
advantages  when  one  wishes  to  calculate  the  "near  field" 
region  or  to  include  shear  wave  effects  in  elastic  media.1112 
It  is  also  often  used  as  a benchmark  for  other  less  exact 
techniques.  We  show  below  the  relationship  between  the 


Rather  than  evaluate  Eq.  (32)  for  each  g for  the 
complete  set  of  ks  [typically  solving  Eq.  (32)  1024  to  8196 
times],  one  can  utilize  a normal  mode  expansion  of  the  form 

g(k,z)  = £ an  (k)un  (z),  (35) 

where  the  quantities  un  are  eigenfucntions  of  the  following 
eigenvalue  problem: 


~ + [K\z)  - kl]  ujtz)  = 0; 

r UJ^L  dz -- 

Jo  p(z) 


where  6^  = 1 for  n = m and  is  zero  for  n ^ m.  The 
eigenfunctions  un  are  zero  at  z = 0,  satisfy  the  local 
boundary  conditions  descriptive  of  the  ocean  bottom 
properties,  and  satisfy  a radiation  condition  for  z -*  oo. 
They  form  an  orthonormal  set  in  a Hilbert  space  with 
weighting  function  p(z),  the  local  density.  Thus,  for 
example,  in  a Pekeris  waveguide,  K is  w/clt2,  where  c1>2  is 
the  speed  of  sound  in  the  water  and  bottom,  respectively, 
and  Eq.  (36)  is  solved  with  the  boundary  conditions  given 
by  Eqs.  (20)  and  (21).  For  the  eigenvalues  and 
eigenfunction,  the  solution  is  Eqs.  (13)  and  (15)  [with  an 
additional  normalization  factor  for  Eq.  (15)  resulting  from 
the  orthonormality  condition  of  Eq.  (36)]. 

The  range  of  discrete  eigenvalues  is  given  by  the 
condition 

min[£(z)]<  kn<  max[A:(z)].  (37) 
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This  is  the  nonisovelocity  generalization  of  Ineq.  (14). 
These  discrete  eigenvalues  correspond  to  discrete  angles 
within  the  critical  angle  cone  in  Fig.  3(a)  such  that  specific 
waves  constructively  interfere.  The  eigenvalues  kn  typically 
have  a small  imaginary  part  a„,  which  serves  as  a modal 
attenuation  coefficient  representative  of  all  the  losses  in  the 
ocean  environment  (see  Ref.  3 for  the  formulation  of 
normal  mode  attenuation  coefficients,  and  other  refer- 
ences). Solving  Eq.  (32)  by  using  the  normal  mode 
expansion  given  by  Eqs.  (35)  and  (36)  permits  an  analytic 
integration  of  Eq.  (3 1 ) yielding  (for  the  source  at  the  origin) 

G{r,z)  = L(z)  £ un(z)un{z)Hl(knr).  (38) 

^ n 

The  asymptotic  form  of  the  Hankel  function  [given  above 
Eq.  (33)]  can  be  used  in  the  above  equation  to  obtain  the 
well-known  normal  mode  representation  of  a cylindrical 
(axis  is  depth)  waveguide: 


Hence,  we  note  that  the  continuous  spectrum  is  the  near 
(vertical)  field  and  the  discrete  spectrum  is  the  [more 
horizontal,  profile-dependent)  far  field  [falling  within  the 
cone  in  Fig.  3(a)]. 

The  advantages  of  the  NM  procedure  are  that: 

• the  solution  is  available  for  all  source  and  receiver 
configurations  once  the  eigenvalue  problem  is 
solved; 

• it  is  easily  extended  to  moderately  range-dependent 
conditions  by  using  the  adiabatic  approximation;  and 

• it  can  be  applied  (with  more  effort)  to  extremely 
range-dependent  environments  using  coupled  mode 
theory. 

However,  it  does  not  include  a full  representation  of  the 
near  field. 


G(r,z) 


W) 

(%vr)m 


exp(-i7r/4) 


X £ TJJi eXP^/  - «/)’ 


(39) 


with  the  modal  attenuation  coefficient  given  by 

„ , » K(z)  | un(z)  | 2 

" K \ o p(z ) 


(40) 


Equation  (39)  is  a far  field  solution  of  the  wave  equation 
and  neglects  the  continuous  spectrum  ( kn  < min[£(z)]  of 
Ineq.  (37)]. 

To  illustrate  the  various  portions  of  the  acoustic  field, 
we  note  that  kn  is  a horizontal  wave  number  so  that  a "ray 
angle"  associated  with  a mode  with  respect  to  the  horizontal 
can  be  taken  to  be  6 = cos1  [kJK(z)].  For  a simple 
isovelocity  waveguide,  the  maximum  sound  speed  is  the 
bottom  sound  speed  corresponding  to  min  [AT(z)].  At  this 
value  of  K(z),  we  have  from  Snell’s  law  6 = 0C,  the  bottom 
critical  angle.  In  effect,  if  we  look  at  a ray  picture  of  the 
modes,  the  continuous  portion  of  the  mode  spectrum 
corresponds  to  rays  with  grazing  angles  greater  than  the 
bottom  critical  angle  of  Fig.  2(b)  and  therefore  outside  the 
cone  of  Fig  3(a).  This  portion  undergoes  severe  loss. 


Adiabatic  and  Coupled  Mode  Theory 

All  of  the  range-independent  normal  mode  "machinery" 
developed  for  environmental  ocean  acoustic  modeling 
applications  can  be  adapted  to  mildly  range-dependent 
conditions  by  using  adiabatic  mode  theory.  The  underlying 
assumption  is  that  individual  propagating  normal  modes 
adapt  (but  do  not  scatter  or  "couple"  into  each  other)  to  the 
local  environment.  The  coefficients  of  the  mode  expansion, 
a in  Eq.  (35),  now  become  mild  functions  of  range,  i.e.,  an 
(k)  -»  an  (k,  r).  This  modifies  Eq.  (39)  as  follows: 


G(r,z) 


ip(z) 

(8*r)1'2 


exp(-ix/4) 


x£ 


F72 


e\p(iTnr  - ocnr). 


(41) 


where  the  range-averaged  eigenvalues  and  attenuation 
coefficients  are 


and  the  kn{r’),  an(r ')  are  obtained  at  each  range  segment 
from  the  eigenvalue  problem  [Eq.  (36)]  evaluated  at  the 
environment  at  that  particular  range  along  the  path.  The 
quantities  un  and  v„  are  the  sets  of  modes  at  the  source  and 
the  field  positions,  respectively. 
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Simply  stated,  the  adiabatic  mode  theory  leads  to  a 
description  of  sound  propagation  such  that  the  acoustic  field 
is  a function  of  the  modal  structure  at  both  the  source  and 
the  receiver  and  some  average  propagation  between  the  two. 
Thus,  for  example,  when  sound  emanates  from  a shallow 
region  where  only  two  discrete  modes  exist  and  propagates 
into  a deeper  region  with  the  same  bottom  (same  critical 
angle),  the  two  modes  from  the  shallow  region  adapt  to  the 
form  of  the  first  two  modes  in  the  deep  region.  However, 
the  deep  region  can  support  many  more  modes.  Intuitively, 
we  therefore  expect  that  the  resulting  two  modes  in  the  deep 
region  will  take  up  a smaller  more  horizontal  part  of  the 
j cone  of  Fig.  3(a)  than  they  take  up  in  the  shallow  region. 
This  means  that  sound  rays  going  from  shallow  to  deep  tend 
to  become  more  horizontal,  which  is  consistent  with  a ray 
picture  of  downslope  propagation. 

Finally,  fully  coupled  mode  theory16  for  range- 
dependent  environments  has  been  developed  but  requires 
additional  intensive  computation.  The  computational  burden 
originates  from  having  to  include  coupling  in  the  backward 
as  well  as  forward  direction.  Under  certain  conditions, 
restricting  the  computation  to  forward  coupling  results  in 
a computationally  tractable  procedure.  For  example, 
forward  coupled-mode  theory  has  been  applied  to  model 
pulse  propagation  for  shallow-water  tomography  studies. 17 

PARABOLIC  EQUATION  MODEL  (PE) 

Standard  PE-Split  Step  Algorithm 


— * * K%J  = 0;  (44) 

dr 2 r dr 

+ 9V  + fj_  + 2 37 
dr2  dz2  [r  7 3rJ  (45) 

+ ?£  + Kin2  t - K&  = 0. 

dr 

Equation  (44)  is  a Bessel  equation,  and  we  take  the 
outgoing  solution,  a Hankel  function  H0l(£or),  in  its 
asymptotic  form  [given  above  Eq.  (33)]  and  substitute  it 
into  Eq.  (45),  together  with  the  "paraxial”  (narrow  angle) 
approximation 

*+  < 2K0  (46) 

dr2  0 dr 

to  obtain  the  parabolic  equation  (in  r), 

H * 2 IK,  * K.V  - W - 0,  (47) 

dz  dr 

where  we  note  that  n is  a function  of  range  and  depth. 
Here,  we  use  the  "split-step"  marching  solution18  to  solve 
the  parabolic  equation.  We  take  n to  be  a constant  for  each 
range  step;  the  error  this  introduces  can  be  made  arbitrarily 
small  by  the  appropriate  numerical  gridding.  The  Fourier 
transform  of  ^ can  then  be  written  as 


The  PE  method  is  a practical  wave-theoretic  range- 
dependent  propagation  model.  In  its  simplest  form,  it  is  a 
farfield  narrow-angle  (—  ±15  deg  with  respect  to  the 
horizontal— adequate  for  many  underwater  propagation 
problems)  approximation  to  the  wave  equation.  Assuming 
azimuthal  symmetry  about  a source,  we  express  the  solution 
of  Eq.  (27)  in  cylindrical  coordinates  in  a source-free  region 
in  the  form 


G(r,z)=  rKr,zyJ(r] ), 


X(r,s)  = _L  [ “ t(r,z)  exp (~isz)dz  (48) 

which  together  with  Eq.  (47)  gives 

-S*x  + 2 iK0  + Klin2  - 1)  X = 0.  (49) 

dr 

The  solution  of  Eq.  (49)  is  simply 


and  we  define  K2  ( r,z ) = K2n2,  n therefore  being  an  "index 
of  refraction"  c0/c,  where  c0  is  a reference  sound  speed. 
Substituting  Eq.  (43)  into  Eq.  (27)  in  a source-free  region 
and  taking  K2  as  the  separation  constant,  7 and  t satisfy  the 
following  two  equations: 


x(r,s ) = x(r0,s) 


x exp 


Kl(n 2 - 1 ) - s2 
2iK0 


(r  - rQ) 


(50) 
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with  specified  initial  condition  at  r0.  The  inverse  transform 
gives  the  field  as  a function  of  depth, 


i{r,z)  = | ” X(r0,s)  exp 


2 


\)Ar 


(51) 


where  zs  is  the  source  depth.  Another  common  method  is  to 
initialize  the  field  with  a normal  mode  representation.  There 
are  also  more  accurate  and  elegant  methods  to  initialize  the 
PE  (refer  to  Chapter  6 of  Ref.  3 for  an  assortment  of 
methods). 

Generalized  or  Higher  Order  PE  Methods 


x exp 


2*o 


exp  (isz)ds , 


where  A r = r - r0.  Introducing  the  symbol  y for  the 
Fourier  transform  operation  from  the  z-domain  [as 
performed  in  Eq.  (48)]  and  y1  as  the  inverse  transform, 
Eq.  (51)  can  be  summarized  by  the  range-stepping 
algorithm, 


Methods  of  solving  the  parabolic  equation,  including 
extensions  to  higher  angle  propagation,'9'22  elastic  media,23 
and  direct  time-domain  solutions  with  nonlinear  effects24 
have  recently  appeared. 

Equation  (47),  with  the  second-order  range  derivative 
that  was  neglected  because  of  Ineq.  (46),  can  be  written  in 
operator  notation  as 


\p(r  + A r,z)  = exp 


iKr 


0 ( n 2 - l)Ar 


x y 


| exp  (- 


2 

iAr 


(52) 


2K° 


s2)  • y m,z)] 


which  is  often  referred  to  as  the  "split  step"  marching 
solution  to  the  PE.  The  Fourier  transforms  are  performed 
using  FFTs.  Equation  (52)  is  the  solution  for  n constant,  but 
the  error  introduced  when  n (profile  or  bathymetry)  varies 
with  range  and  depth  can  be  made  arbitrarily  small  by 
increasing  the  transform  size  and  decreasing  the  range-step 
size. 


Finally,  since  the  PE  is  an  initial  value  problem,  we 
must  have  a starting  field.  Typically  we  can  approximate  a 
point  source  as  a Gaussian  shape  since  we  are  only 
concerned  with  propagation  angles  confined  to  the  narrow 
cone  discussed  in  this  paper:3 


mz)  = 


(53) 


When  the  source  is  near  the  surface,  it  is  helpful  to  use  an 
antisymmetric  combination  of  the  source  and  its  image  so 
that  the  surface  boundary  condition  is  satisfied, 


[J P 2 + 2iK0P  + Ko(Q 2 - 1)]  </-  = 0,  <55) 


where 


P = 


+ ±$_ 
+ K20dz2' 


(56) 


Factoring  Eq.  (55)  assuming  weak  range-dependence  and 
retaining  only  the  factor  associated  with  outgoing  prop- 
agation yields  a one-way  equation 

P'P  = iK0(Q  - l)i,  (57) 

which  is  a generalization  of  the  parabolic  equation  beyond 
the  narrow  angle  approximation  associated  with  Ineq.  (46). 

If  we  define  Q = yjl  + q and  expand  Q in  a Taylor  series 
as  a function  of  q , the  standard  PE  method  is  recovered 
with  Q * 1 + .5^.  The  wide-angle  PE  to  arbitrary 
accuracy  in  angle,  phase,  etc.,  can  be  obtained  from  a Pade 
series  representation  of  the  Q operator,19 

(2  = v/w  =l+£-J%^  + 0(«?2'-'),  (58) 

ft  1 


where  m is  the  number  of  terms  in  the  Pade  expansion  and 


<A(0,z)  = iK0,z  - z)  - M0,Z  + z),  (54) 


a. 


2 

2m+l 


sin2 


2m+\ 


fi.  = cos2 

J,m 


Jjj 

2m+l 


(59) 
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The  solution  of  Eq.  (57)  using  Eqs.  (58)-(59)  has  been 
implemented  using  finite-difference  techniques  for  fluid  and 
elastic  media.19,23 

A numerical  method  that  allows  one  to  solve  the  PE 
with  large  range  steps  (although  the  environmental  range- 
dependence  must  still  be  adequately  sampled)  has  been 
recently  developed  and  is  referred  to  as  the  split-step 
Pade.20,21  In  this  technique,  Eq.  (57)  is  integrated 
analytically, 

P(r + A r)  = exp  [zX0(  - 1 + Jl+q  ’ )]P(r) , (^0) 

and  the  the  Pade  approximation  is  then  applied, 

W SJ  Q 

exp[i*o(-l  + y[TTJ)]  = £ , (61) 

U 1 + h«a 

where  the  coefficients  in  the  Pade  series  are  discussed  in 
Ref.  21.  The  split-step  Pade  formula  is  obtained  by 
substituting  Eq.  (61)  into  Eq.  (60), 


P(r+Ar)  = P(r)  + £ (a.m)(l  +bjmq)-'qP(r).  (62) 

Jm  1 

! 

Range  steps  of  more  than  an  order  of  magnitude  greater 
than  that  used  in  the  previous  Pade  implementation  are  now 
permissible.  Furthermore,  the  structure  of  the  solution 
permits  the  efficient  use  of  parallel  processors  to  speed  up 
execution  of  the  algorithm.  As  far  as  the  starting  fields  are 
concerned,  the  same  issues  apply  as  discussed  in  the 
previous  subsection.  However,  the  "self-starter"  algorithm22 
is  particularly  useful. 

Transmission  Loss  (TL) 

The  models  discussed  in  this  section  compute  complex 
quantities  such  as  pressure.  However,  in  underwater 
acoustics  it  is  convenient  to  express  the  field  in  terms  of 
transmission  loss,  which  is  defined  as 

i 

TLd-.r,)  = -20  log10  l^yl.  <63> 
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where  p( r,  r,)  is  the  acoustic  pressure  at  point  r from  a 
simple  point  source  at  point  rJ7  and  p0(r 0)  is  the  pressure 
produced  at  a distance  of  1 m from  the  same  source  in  an 
infinite  homogeneous  medium  with  density  p(rJ).  Thus,  we 
have  that  p0( r0)  = (47rr)1  exp (iKr)  with  r = 1 . This  formula 
is  also  referred  to  as  the  coherent  TL  when  applying  it  to 
normal  mode  models  in  that  TL  contains  the  variations 
resulting  from  interfering  modes.  The  incoherent  TL  is 
obtained  from  a normal  mode  model  by  taking  for  the 
numerator  the  square  root  of  the  sum  of  the  magnitudes 
squared  of  the  individual  modal  terms  in  Eqs.  (39)  or  (41). 

Propagation  of  Pulses 

Presendy,  the  most  practical  method  to  compute  the 
shape  of  a pulse  as  it  propagates  in  shallow  water  is  by 
Fourier  synthesis.  That  is,  a transmitted  pulse  with  a source 
spectrum /(w)  will  have  a pulse  shape  at  (r,z,0  given  by 

P(r,z,t)  = | “ j{u>)p(j,z,u>)  e'ku  da,  (64) 

where  p(r9z;o))  is  the  acoustic  field  that  can  be  computed 
from  any  of  the  wave  models  discussed  in  this  section.  In 
particular,  a normal  mode  over  moderate  bandwidths  can  be 
used  very  efficiendy  by  interpolating  modal  wavenumbers 
across  frequency  but  assuming  that  the  mode  shapes  do  not 
significantly  change.  This  interpolation  has  even  been 
applied  to  coupled-mode  computations.17  Time  series  can 
also  be  generated  efficiendy  using  ray  theory.25 

SOME  RESULTS 

Range-Independent  TL 

It  is  well  established  that  all  the  models  should  be  in 
agreement  for  range-independent  environments.  In  over- 
lapping regions  of  validity  (far  field),  normal  mode  and 
spectral  (FFP)  solutions  are  identical.  It  is  only  more 
recendy  the  PE  results  have  converged  to  give  the  same 
solution  as  a benchmark  FFP  solution.  Figure  5 shows  some 
results3  for  a 100-m  Pekeris  waveguide  with  source  and 
receiver  at  a depth  of  99.5  m (for  maximum  interference 
structure).  The  water  and  bottom  speeds  are  1500  and  1590 
m/s,  respectively;  the  density  ratio  between  water  and 
bottom  is  1.2;  and  the  bottom  attenuation  is  0.5  dB/X.  The 
frequency  is  250  Hz.  The  standard  PE  without  high  angle 
accuracy  fails  to  match  the  benchmark  FFP  results  whereas 
the  higher  angle  PE  does. 
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Fig.  5 — Comparison  of  transmission  loss  results  for  narrow 
and  wide-angle  PEs  with  FFP  reference  solution 


Bottom  Shear  Wave  Effects 

For  ocean  bottoms  with  at  least  some  rigidity,  the 
effect  of  the  existence  of  shear  provides  the  sound  field  with 
additional  degrees  of  freedom  to  interact  with  the  bottom. 
The  mechanism  is  coupling  into  shear  waves  and/or 
interface  waves  (for  a source  acoustically  close  to  the 
bottom  interface).  For  bottoms  with  shear  wave  speeds  less 
than  the  speed  of  sound  in  the  water  column,  the  shear 
waves  carry  energy  away  from  the  water  column.  The 
major  effect  is  increased  attenuation  at  low  frequencies  over 
nonshear- supporting  ocean  bottoms.  For  the  case  when  the 
shear  speed  is  higher  than  the  water  column  sound  speed, 
the  ratio  of  water  column  sound  speed  to  shear  speed 
determines  the  critical  angle.  If  a source  is  close  (of  the 
order  of  a wavelength)  to  the  bottom,  it  will  couple  into  an 
interface  wave  (Scholte  wave  at  fluid-solid  interface  or 
Stoneley  wave  and  solid-solid  interface).  These  wave  can  be 
detected  by  hydrophones  near  or  on  the  bottom  or 
seismometers  on  the  bottom.  All  these  effects  are 
computable  using  an  FFP  model.12 

Broadband  Modeling 

Because  modes  travel  at  different  phase  and  group 
speeds  as  a function  of  frequency,  a pulse  will  disperse  as 
it  propagates.  The  phase  speed  is  the  horizontal  propagation 
speed  ( un  = o)!kn)  of  a wavefront  [see,  for  example  Fig. 
3(b)]  corresponding  to  a mode,  whereas  the  group  speed  (vn 
- (< dkjdo) )'1)  of  a mode  represents  the  speed  at  which  an 
energy  packet  (which  has  a finite  bandwidth)  propagates 
down  the  waveguide.  Figure  6 shows  plots  of  modal  group 
speeds  and  a dispersing  pulse  in  a Pekeris  waveguide.3  The 
sound  speeds  in  the  water  column  and  bottom  are  1500  and 


(a)  Dispersion  characteristics  of  the  waveguide 


(c)  Received  signal  at  a distance  of  30  km  from  the  source. 


Fig.  6 — 50  Hz  bandwidth  pulse  propagation  in  a Pekeris  wave- 
guide.3 The  initial  pulse  length  of  0.08  s becomes  strongly  dispersed, 
with  a final  length  of  around  0.7  s. 
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1600  m/s,  respectively,  The  density  ratio  is  1.5,  and  the 
bottom  attenuation  is  0.2  dB/X.  Here,  we  see  that  the  initial 
pulse  disperses  into  two  distinct  pulses  (mode  1 and  2)  as  a 
result  of  the  different  modal  group  speeds.  At  the  trailing 
tail  is  a weak  third  mode  arrival. 

I 

We  also  should  note  that,  in  general,  a good 
j approximation  to  modeling  TL  over  some  bandwidth,  say, 
a one- third  octave  band,  is  simply  to  compute  the  incoherent 
normal  mode  transmission  loss  at  the  center  frequency. 

I 

Optimum  Frequency 

The  tradeoff  between  high-frequency  attenuation  and 
scattering  and  low-frequency  bottom  loss  leads  to  an 
optimum  frequency  of  propagation  in  a waveguide.  The 
optimum  frequency  in  a shallow-water  environment  depends 
to  a large  extent  on  the  width  of  the  effective  ducted 
j propagation  as  well  as  the  properties  of  the  water  column 
and  bottom  type.  For  a winter  profile,  the  duct  corresponds 
to  the  water  depth.  For  a summer  profile  with  source  and 
receiver  below  the  thermocline,  the  effective  duct  is  the 
water  column  below  the  thermocline.  Figure  7 is  the  results 
of  an  experimental  and  modeling  study  including  a contour 
of  third-octave  TL  vs  frequency  and  range.26  The  optimum 
j frequency  is  seen  to  be  about  in  the  region  of  200-400  Hz 
from  the  contours.  For  example,  at  a range  of  about  60  km, 
the  loss  increases  for  frequencies  above  and  below  the 
optimum  frequency  of  propagation.  Agreement  between 
experiment  and  theory  was  obtained  by  using  the  model/data 
comparison  to  fine-tune  the  bottom  parameters.  Of  course, 

] one  set  of  bottom  parameters  must  satisfy  the  data  at  all 
j frequencies. 

j 

Range-Dependent  TL 

The  PE  is  the  most  robust  method  for  computing  mid- 
to  low-frequency  TL  in  a range-dependent  shallow-water 
environment.  The  April  1990  issue  of  the  Journal  of  the 
Acoustical  Society  of  America  has  a set  of  papers 
benchmarking  various  range-dependent  models.27  A range- 
dependent  phenomenon  of  growing  interest  is  the  existence 
of  internal  wave  solitons.  These  have  been  studied 
deterministically  to  demonstrate  their  acoustic  influence. 
Using  the  parabolic  equation  and  projecting  the  results  into 
mode  space,  it  has  been  shown  that  these  solitons  can 
explain  measured  data  in  such  an  environment.28 

Three-Dimensional  (3D)  Propagation 
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approximation.  There  are  essentially  two  methods  to  treat 
a three-dimens ionally  varying  environment.  The  most 
straightforward  and  still  not  very  practical  method  is  to 
formally  solve  the  3D  wave  equation  for  a 3D  environment. 

This  has  been  done  to  some  extent  using  ray  methods29  but 
is  still  rather  numerically  intractable  in  the  wave  theoretic 
regime.  An  approximation  is  to  compute  2D  TL  along 
radials  and  outward  from  a source.30  Figure  8 shows  an 
example  of  this  for  a bathymetrically  variable  shallow- water  \ 

region.31  TL  loss  computations  along  horizontal  refracting 
paths  have  also  been  made  using  a new  model  referred  to  as 
the  Adiabatic  Mode  PE.32  This  latter  reference  includes 
very-low-frequency  deep-water  computations  that  essentially 
scale  to  approximately  the  shallow-water  acoustic  regime. 

SUMMARY 

An  assortment  of  models  have  been  developed  to 
compute  the  acoustic  field  in  shallow-water  environments  in 
both  the  frequency  and  time  domains.  This  paper  has  not 
considered  fluctuations  and  other  forms  of  statistical  spatio- 
temporal  variability  of  the  environment/acoustics.  The 
relationship  between  these  models  is  well  understood,  and 
their  accuracy  as  related  to  solving  the  wave  equation  is 
also  understood.  However,  these  models  require 
environmental  inputs— sound  speeds  and  geophysical 
properties  of  the  bottom  as  a function  of  position— which 
turn  out  to  be  the  limiting  factor  in  how  well  they  model 
real -world  scenarios. 
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Fig.  7 — Optimum  frequency,  (a)  Summer  shallow-water  environment;  (b)  Frequency  vs  depth  propagation  loss  contours; 
experiment  and  model  results;  (c)  Third  octave  propagation  loss  curves  (which  are  horizontal  cuts  of  contours 
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Fig.  8 — Transmission  loss  map  constructed  from  adiabatic  mode  model 
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Appendix 

UNITS 


} 

The  decibel  (dB)  is  the  dominant  unit  in  underwater 
acoustics  and  denotes  a ratio  of  intensities  (not  pressures) 
expressed  in  terms  of  a logarithmic  (base  10)  scale.  Two 
intensities,  I{  and  /2,  have  a ratio,  /j//2,  in  decibels  of  10  log 
I{/I2  dB.  Absolute  intensities  can  therefore  be  expressed  by 
using  a reference  intensity.  The  presently  accepted  reference 
intensity  is  a micropascal  (/xPa):  the  intensity  of  a plane 
wave  having  an  rms  pressure  equal  to  10'5  dynes  per  square 
! centimeter.  Therefore,  taking  1 jiPa  as  /2,  a sound  wave 
having  an  intensity,  of,  say,  one  million  times  that  of  a 
plane  wave  of  rms  pressure  1 /iPa  has  a level  of  10  log 
(106/1)  m 60  dB  re  1 /xPa.  Pressure  (p ) ratios  are  expressed 
in  dB  re  1 /uPa  by  taking  20  log  pxlp2,  where  it  is 
understood  that  the  reference  originates  from  the  intensity 
of  a plane  wave  of  pressure  equal  to  1 /iPa. 
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